#### Table S16: Multilevel regression (only 2023) ####

# Libraries
# library(here)
# library(rio)
# library(tidyverse)
# library(lubridate)
# library(stargazer)
# library(lme4)

# data_pnas = import(here("Data","data_pnas.rds"))

p_overturn_mod1 = lmer(norm_loyaltyre ~ voted_overturn + (1 | district) + (1 | uid), 
                       data = data_pnas |> filter(date >= mdy("01/03/2023")), weights = weight)
p_overturn_mod2 = lmer(norm_loyaltyre ~ voted_overturn*pid + (1 | district) + (1 | uid), 
                       data = data_pnas |> filter(date >= mdy("01/03/2023")), weights = weight)
p_denial_538_mod1 = lmer(norm_loyaltyre ~ denied_538 + (1 | district) + (1 | uid), data = 
                           data_pnas |> filter(date >= mdy("01/03/2023")), weights = weight)
p_denial_538_mod2 = lmer(norm_loyaltyre ~ denied_538*pid + (1 | district) + (1 | uid),  
                         data = data_pnas |> filter(date >= mdy("01/03/2023")), weights = weight)

stargazer(p_overturn_mod1, p_denial_538_mod1, p_overturn_mod2, p_denial_538_mod2,
          title = "Multilevel Regression Results - 2023 Only",
          dep.var.labels = "Loyalty",
          covariate.labels = c("Voted to Overturn",
                               "Election Denial",
                               "Independent",
                               "Republican",
                               "Voted to Overturn:Independent",
                               "Voted to Overturn:Republican",
                               "Election Denial:Independent",
                               "Election Denial:Republican"),
          notes = "Estimated with survey weights and two-sided tests",
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat = c("aic","bic"),
          type = "latex",
          out = here("Tables","Supplementary","table_s16.tex"),
          header = F)